*FIG 9: EVENT STUDIES OF TAX REVENUE

/* following packages need to be installed:
ssc install winsor2
ssc install xtvent
ssc install coefplot
*/



version 16.1

clear all
cap log close 
cap clear matrix
set more off
graph drop _all
set scheme mygraphs

cd "$mypathRR"

* SET UP CANTONAL TAX REVENUES AND POPULATION DATA
*use "Datasets/rev_ktn+gdn_1990-2016.dta", clear


use "Datasets/revenuedata_ktn+gde.dta", clear

merge 1:1 year canton using "Datasets/pop_1971-2017.dta"
drop if year < 1990
//drop if year == 2015 // because of huge one-time-effect in Sarnen

* Label & rescale variables *
label var year "Year"

* attach canton value labels
run "Resources/labels-cantons.do"
label values canton cant


// correct the outliers by winsorizing
	* winsorize observations with a high revenues (it's all about income in 2015 - only winsorizing at 96th percentile gets rid of that outlier)
	foreach v in  rev_corp_total revCH_corp_total revCentralOW_corp_total revCentral_corp_total rev_total rev_inctax rev_weatax rev_sourcetax rev4008 rev4009 revCH_total revCH_inctax revCH_weatax revCH_sourcetax revCH4008 revCH4009 revCentral_total revCentral_inctax revCentral_weatax revCentral_sourcetax revCentral4008 revCentral4009 revCentralOW_total revCentralOW_inctax revCentralOW_weatax revCentralOW_sourcetax revCentralOW4008 revCentralOW4009 {
	winsor2 `v' , cuts(0 99) suffix(_w) by(cant)
	rename `v' `v'_raw
	rename `v'_w `v'
	}


* generate p.c. variables
foreach ending in _total _inctax _weatax _sourcetax 4008 4009 _corp_total {
  gen pcrev`ending'=rev`ending'/totpop
  gen pcrevCH`ending'=revCH`ending'/totpopCH
  gen pcrevCentral`ending'=revCentral`ending'/popCentral
  gen pcrevCentralOW`ending'=revCentralOW`ending'/popCentralOW
}


// label p.c. variables
foreach v in  rev_corp_total revCH_corp_total revCentralOW_corp_total revCentral_corp_total rev_total rev_inctax rev_weatax rev_sourcetax rev4008 rev4009 revCH_total revCH_inctax revCH_weatax revCH_sourcetax revCH4008 revCH4009 revCentral_total revCentral_inctax revCentral_weatax revCentral_sourcetax revCentral4008 revCentral4009 revCentralOW_total revCentralOW_inctax revCentralOW_weatax revCentralOW_sourcetax revCentralOW4008 revCentralOW4009 {
	local u : variable label `v'
	local l= "`u' p.c." 
	label var pc`v' "`l'"
  }
  
  
* gen log(revenue) in Mio CHF
foreach var in rev_corp_total revCH_corp_total revCentralOW_corp_total revCentral_corp_total  rev_total revCH_total revCentralOW_total revCentral_total rev_inctax revCH_inctax revCentralOW_inctax revCentral_inctax rev_weatax revCH_weatax revCentralOW_weatax revCentral_weatax rev_sourcetax revCH_sourcetax revCentralOW_sourcetax revCentral_sourcetax rev4008 revCH4008 revCentralOW4008 revCentral4008 rev4009 revCH4009 revCentralOW4009 revCentral4009 {
gen ln`var' = .
replace ln`var' = log(`var'/1000)

// label logged revenue variables
	local u : variable label `var'
	local l= "`u' (log)" 
	label var ln`var' "`l'"
  }
  
  rename canton cant
*--------------------------------------------------------------------------------------
* PREPARE DATA FOR EVENT STUDY

cd "$mypathRR/Results"

* Gen canton treatment dummies
gen treated = (cant == 6)
  label var treated "Treated"
  
gen period = (year > 2005 & year!=.)
  label var period "Period $ t>2005$"

* Gen DiD Interaction term
gen interaction = treated*period
  label var interaction "DiD"


* Gen interaction terms for pre-treatment periods
forval n=1/16 {
local i= 2006 -`n'
gen pre`n'=0
replace pre`n'= 1 if treated == 1 & year == `i'
label var pre`n' "`i' (treatment lag `n')"
}
replace pre1 = 0 // make 2005 the reference year


* Gen interaction terms for post-treatment periods
forval n=1/11 {
local i = 2005 + `n'
gen post`n' = 0
replace post`n' = 1 if treated == 1 & year == `i'
label var post`n' "`i' (treatment lead `n')"
}


* Gen canton-specific time trends
	tab cant, gen (cant_d)
	// adjust names for correct canton id
	local i = 0
	forval n = 0/26 {
		local i = `i'+1
		rename cant_d`i' cant_d`n'
	}
	
	foreach var of varlist cant_d* {
	gen trend_`var' = year*`var'
	label var trend_`var' "Canton specific time trend"
	}
		
	// overall trend
	gen trend = year


* * * * define the correct weight variable for each outcome * * * * 

* get number of taxpayers in pre-treatment period
	frame put totpop totpopCH popCentral popCentralOW cant year , into(weights1)
	frame change weights1
	
		keep if year == 2005
		drop year
		foreach var in totpop totpopCH popCentral popCentralOW {
			rename `var' `var'_2005
		}
		
	frame change default

* link back to main frame
	frlink m:1 cant, frame(weights1) gen(link_weights1)
	frget totpop_2005 totpopCH_2005 popCentral_2005 popCentralOW_2005, from(link_weights1)

* replace with actual population in pre-treatment period
	foreach var in totpop totpopCH popCentral popCentralOW {
	replace `var'_2005 = `var' if year < 2005
	}

* generate a weigth variable for each outcome

foreach ending in _total _inctax _weatax _sourcetax 4008 4009 _corp_total {
  gen weight_pcrev`ending'			= totpop_2005
  gen weight_pcrevCH`ending'		= totpopCH_2005
  gen weight_pcrevCentral`ending'	= popCentral_2005
  gen weight_pcrevCentralOW`ending'	= popCentralOW_2005
}

* the total revenue regressions should not be weighted, hence weight == 1
foreach var in rev_corp_total revCH_corp_total revCentralOW_corp_total revCentral_corp_total  rev_total revCH_total revCentralOW_total revCentral_total rev_inctax revCH_inctax revCentralOW_inctax revCentral_inctax rev_weatax revCH_weatax revCentralOW_weatax revCentral_weatax rev_sourcetax revCH_sourcetax revCentralOW_sourcetax revCentral_sourcetax rev4008 revCH4008 revCentralOW4008 revCentral4008 rev4009 revCH4009 revCentralOW4009 revCentral4009 {
	 gen weight_ln`var'	= 1	 
	 gen weight_`var'	= 1

  }
  
  
  
* * *   * * *   * * *   * * *   * * *   * * *
*  RESIDUALIZED  TIME TREND ESTIMATION. * PROBABLY NOT NEEDED FOR PUBLICATION
* * *   * * *   * * *   * * *   * * *   * * *
global outcomes "pcrev_corp_total pcrev_total pcrev_inctax  pcrev_weatax lnrev_corp_total lnrev_total lnrev_inctax  lnrev_weatax rev_total rev_inctax  rev_weatax"

* Estimate a time trend for each  outcome	 	
* estimate & predict canton trend
	preserve
		foreach outcome in $outcomes  {
			reg `outcome' trend_* cant_d* if year < 2006 & year > 1994 [aweight = weight_`outcome']
			predict tr`outcome', xb
			label var tr`outcome' "predicted cantonal time trend in `outcome'"
		}
		
		keep  cant year tr* 
		drop trend* treated
		bys cant year: keep if _n==1
	save "est_canton_trends.dta", replace
	restore 	 
	drop cant_d*

* match trend back to original data
	merge m:1 year cant using "est_canton_trends.dta", nogen keep(1 3)

	rm  "est_canton_trends.dta"

* generate residualized outcomes
foreach outcome in $outcomes {
	reg `outcome' tr`outcome'
	predict resid`outcome', residual
	label var resid`outcome' "residualized `outcome' (after taking out canton trend)"
	
	gen detrend`outcome' = `outcome' - tr`outcome'
	label var detrend`outcome' "detrended `outcome' (after taking out canton trend)"
}


* check residualized outcomes
/*
scatter residpcrev_inctax detrendpcrev_inctax if cant == 6 || lfit residpcrev_inctax detrendpcrev_inctax if cant == 6, xlab(-0.5(0.25)2, grid) ylab(-0.5(0.25)2, grid)
scatter residpcrev_weatax detrendpcrev_weatax if cant == 6 || lfit residpcrev_weatax detrendpcrev_weatax if cant == 6, xlab(-0.25(0.1).25, grid) ylab(-0.25(0.1).25, grid)

scatter pcrev_weatax year if cant == 6 || lfit pcrev_weatax year if cant == 6 // positve trend in wealth
scatter residpcrev_weatax year if cant == 6 || lfit residpcrev_weatax year if cant == 6 // negatvie trend in residualized wealth
scatter detrendpcrev_weatax year if cant == 6 || lfit detrendpcrev_weatax year if cant == 6 // negatvie trend in detrended wealth
scatter detrendpcrev_weatax year if cant == 6 || lfit detrendpcrev_weatax year if cant == 6 & year < 2006 || lfit detrendpcrev_weatax year if cant == 6 & year > 2005 
	// almost flat trend in detrended wealth pre reform

scatter pcrev_inctax year if cant == 6 || lfit pcrev_inctax year if cant == 6 // positve trend in income
scatter residpcrev_inctax year if cant == 6 || lfit residpcrev_inctax year if cant == 6 // only sligthly positive trend in income
scatter residpcrev_inctax year if cant == 6 || lfit residpcrev_inctax year if cant == 6 & year != 2015 // positive slope driven by 2015

*/
*--------------------------------------------------------------------------------------
*--------------------------------------------------------------------------------------
*--------------------------------------------------------------------------------------
* EVENT STUDY GRAPHS with xtevent

cd "$mypathRR/Results/"


xtset cant year

gen policyvar = (cant == 6 & year>2005)

global outcomespaper "pcrev_weatax pcrev_inctax pcrev_total"

// weighted regression, cluster canton-year and gemeinde, 2 models in 1 plot
foreach outcome in $outcomespaper {
	
	// save pre-reform levels to compute percentage changes
	sum `outcome' if cant == 6 & year == 2005
	global `outcome'2005 = `r(max)'
	
	
	// diff-in-diff 1: residualized outcome
	reghdfe resid`outcome' interaction treated period if cant > 0 & cant <27 [aweight = weight_`outcome'] , vce(cluster cant year) absorb(year cant)
	global ES_b_resid`outcome'=_b[interaction]
	global ES_se_resid`outcome'=_se[interaction]
	// percentage change of coefficients
	lincom interaction/${`outcome'2005}*100
	global ES_b_resid`outcome'_pct=`r(estimate)'
	global ES_se_resid`outcome'_pct=`r(se)'
	
	* xtevent model 1: pre-trend correction
	xtevent `outcome'  [aweight = weight_`outcome'] if cant > 0 & cant <27, policyvar(policyvar) 	///
	window(-5 9) trend(-5) reghdfe										///
	vce(cluster cant year)  
	est sto `outcome'_pretrend_II
	
	xteventplot, name(comp_`outcome'II, replace) nosupt minus1label 		///
	ytitle("`: variable label `outcome''") 										///
	ylabel(`myla', labsize(small) grid) 										///
	xlab( , labsize(small) ang(45)) 				///
	xline(-0.2, lcolor(red)) xline(1.8, lcolor(green)) ///
	ciplotopts(recast(rarea) color(maroon*0.2) color(maroon%10))


	// diff-in-diff 2: no pre-trend correction
	reghdfe `outcome' interaction treated period if cant > 0 & cant <27 [aweight = weight_`outcome'] , vce(cluster cant year) absorb(year cant)
	global ES_b_`outcome'=_b[interaction]
	global ES_se_`outcome'=_se[interaction]
	// percentage change of coefficients
	lincom interaction/${`outcome'2005}*100
	global ES_b_`outcome'_pct=`r(estimate)'
	global ES_se_`outcome'_pct=`r(se)'

	* xtevent model 2: no pre-trend correction
	xtevent `outcome'  [aweight = weight_`outcome'] if cant > 0 & cant <27, policyvar(policyvar) 	///
	window(-5 9) reghdfe  													///
	vce(cluster cant year)  
	est sto `outcome'

	xteventplot, name(comp_`outcome', replace) nosupt minus1label  			///
	ytitle("`: variable label `outcome''") 										///
	ylabel(`myla', labsize(small) grid) 										///
	xlab(, labsize(small) ang(45)) 				///
	xline(-0.2, lcolor(red)) xline(1.8, lcolor(green)) ///
	ciplotopts(recast(rarea) color(maroon*0.2) color(maroon%10))
}

	
* 1. ESTIMATES WITHOUT PRE-TREND CORRECTION
/* 	extract estimation matrix / add _k_eq_m5 to _k_eq_m1 to matrix
	estimates stored in e(b) OR e(delta)
	covariance matrix is e(V) OR e(Vdelta)	
*/
	 foreach OUT in weatax inctax total { 
		estimates restore pcrev_`OUT'
		
		* manually calculate confidence intervals from stored results
		loc df = e(df)
		di `df'
		tempname b V
		mat `b' = e(delta)
		mat `V' = e(Vdelta)
		tempvar se	
		
		* standard erors
	    mata st_matrix("`se'", sqrt( diagonal(st_matrix("`V'"))) ')	

		mat list `b'
		mat list `se'
		
		* confidence intervals
		forval n = 1/16 {
		loc ll`n' = `b'[1,`n'] - invttail(`df',0.5*0.05) * `se'[1,`n'] 
		loc ul`n' = `b'[1,`n'] + invttail(`df',0.5*0.05) * `se'[1,`n'] 
		}

		* create a new matrix to put the results
		matrix B`OUT' = J(3,17,.)
        matrix rownames B`OUT' = est ll95 ul95
        matrix colnames B`OUT' = _k_eq_m6 _k_eq_m5 _k_eq_m4 _k_eq_m3 _k_eq_m2 _k_eq_m1 _k_eq_p0 _k_eq_p1 _k_eq_p2 _k_eq_p3 _k_eq_p4 _k_eq_p5 _k_eq_p6 _k_eq_p7 _k_eq_p8 _k_eq_p9 _k_eq_p10
		mat list B`OUT'
		
		* put the results from xtevent into the new matrix
		forval n = 1/5 {
		loc i = `n'
		matrix B`OUT'[1,`n'] = `b'[1,`i'] \ `ll`i'' \ `ul`i''
		}
	 	matrix B`OUT'[1,6] = 0 \ 0 \ 0
		forval n = 7/17 {
		loc i = `n'-1
		matrix B`OUT'[1,`n'] = `b'[1,`i'] \ `ll`i'' \ `ul`i''
		}
		
		* compare results stored in new matrix with output from xtevent
		mat list B`OUT'
		mat list `se'
		estimates replay pcrev_`OUT'
	 }
		
		coefplot (matrix(Bweatax), ci((2 3)) label(" ") omitted color(maroon*1) ciopts(recast(rarea) color(maroon*0.2) color(maroon%10)) ///
			transform(* = @/$pcrev_weatax2005*100)) ///
				(matrix(Binctax), ci((2 3)) label(" ") omitted color(navy*0.6) ciopts(recast(rarea) color(navy*0.2) color(navy%10)) ///
			transform(* = @/$pcrev_inctax2005*100)) ///
			(matrix(Btotal), ci((2 3)) label(" ") omitted color(midblue*1.1) ciopts(recast(rarea) color(midblue*0.3) color(midblue%25)) ///
			transform(* = @/$pcrev_total2005*100)) ///
					,  ///
					coeflabels(  _k_eq_m6 = "-2000" 			///
					_k_eq_m5  = "2001" _k_eq_m4 = "2002" _k_eq_m3 = "2003" _k_eq_m2 = "2004" _k_eq_m1 = "2005" 	///
					_k_eq_p0 = "2006" _k_eq_p1 = "2007" _k_eq_p2 = "2008" _k_eq_p3 = "2009" 	///
					_k_eq_p4 = "2010" _k_eq_p5 = "2011" _k_eq_p6 = "2012" _k_eq_p7 = "2013" 	///
					_k_eq_p8 = "2014" _k_eq_p9 = "2015"  _k_eq_p10 = "2016+") 					///
					vertical levels(95) ciopts(recast(rarea)) connect(direct) nooffsets ///
					yline(0, lcolor(black))  xline(6.4, lcolor(red)) xline(8.7, lcolor(green)) ///
					xlabel(, labsize(small) ang(45)) ylabel(-80(20)60, labsize(small) grid) ///
					ytitle("Impact on revenue per capita (in %)") name(xtevent_NOtrend_comb_percentage, replace) ///
					legend(order(6 "Total personal tax"	- "DiD-estimate: {&beta} = `:di %3.2f ${ES_b_pcrev_total_pct}' (`:di %3.2f ${ES_se_pcrev_total_pct}')" ///
								4 "Income tax" - "DiD-estimate: {&beta} = `:di %3.2f ${ES_b_pcrev_inctax_pct}' (`:di %3.2f ${ES_se_pcrev_inctax_pct}')"  ///
								2 "Wealth tax" -  "DiD-estimate: {&beta} = `:di %3.2f ${ES_b_pcrev_weatax_pct}' (`:di %3.2f ${ES_se_pcrev_weatax_pct}')" ) ///
								row(3) ring(0) position(11) bmargin(small) size(small)) 
								
	
* 2. ESTIMATES WITH PRE-TREND CORRECTION
/* 	extract estimation matrix / add _k_eq_m5 to _k_eq_m1 to matrix
	estimates stored in e(b) OR e(delta)
	covariance matrix is e(V) OR e(Vdelta)	
*/
	 foreach OUT in weatax inctax total { 
		estimates restore pcrev_`OUT'_pretrend_II
		
		* manually calculate confidence intervals from stored results
		loc df = e(df)
		di `df'
		tempname b V
		mat `b' = e(delta)
		mat `V' = e(Vdelta)
		tempvar se	
		
		* standard erors
	    mata st_matrix("`se'", sqrt( diagonal(st_matrix("`V'"))) ')	

		mat list `b'
		mat list `se'
		
		* confidence intervals
		forval n = 1/12 {
		loc ll`n' = `b'[1,`n'] - invttail(`df',0.5*0.05) * `se'[1,`n'] 
		loc ul`n' = `b'[1,`n'] + invttail(`df',0.5*0.05) * `se'[1,`n'] 
		}

		* create a new matrix to put the results
		matrix C`OUT' = J(3,17,.)
        matrix rownames C`OUT' = est ll95 ul95
        matrix colnames C`OUT' = _k_eq_m6 _k_eq_m5 _k_eq_m4 _k_eq_m3 _k_eq_m2 _k_eq_m1 _k_eq_p0 _k_eq_p1 _k_eq_p2 _k_eq_p3 _k_eq_p4 _k_eq_p5 _k_eq_p6 _k_eq_p7 _k_eq_p8 _k_eq_p9 _k_eq_p10
		mat list C`OUT'
		
		* put the results from xtevent into the new matrix
	 	matrix C`OUT'[1,1] = `b'[1,1] \ `ll1' \ `ul1'
	 	matrix C`OUT'[1,2] = 0 \ 0 \ 0
	 	matrix C`OUT'[1,3] = 0 \ 0 \ 0
	 	matrix C`OUT'[1,4] = 0 \ 0 \ 0
	 	matrix C`OUT'[1,5] = 0 \ 0 \ 0
	 	matrix C`OUT'[1,6] = 0 \ 0 \ 0
		forval n = 7/17 {
		loc i = `n'-5
		matrix C`OUT'[1,`n'] = `b'[1,`i'] \ `ll`i'' \ `ul`i''
		}
		
		* compare results stored in new matrix with output from xtevent
		mat list C`OUT'
		mat list `se'
		estimates replay pcrev_`OUT'_pretrend_II
	 }
		

		coefplot (matrix(Cweatax), ci((2 3)) label(" ") omitted color(maroon*1) ciopts(recast(rarea) color(maroon*0.2) color(maroon%10)) ///
			transform(* = @/$pcrev_weatax2005*100)) ///
				(matrix(Cinctax), ci((2 3)) label(" ") omitted color(navy*0.6) ciopts(recast(rarea) color(navy*0.2) color(navy%10)) ///
			transform(* = @/$pcrev_inctax2005*100)) ///
			(matrix(Ctotal), ci((2 3)) label(" ") omitted color(midblue*1.1) ciopts(recast(rarea) color(midblue*0.3) color(midblue%25)) ///
			transform(* = @/$pcrev_total2005*100)) ///
					,  ///
					coeflabels(  _k_eq_m6 = "-2000" 			///
					_k_eq_m5  = "2001" _k_eq_m4 = "2002" _k_eq_m3 = "2003" _k_eq_m2 = "2004" _k_eq_m1 = "2005" 	///
					_k_eq_p0 = "2006" _k_eq_p1 = "2007" _k_eq_p2 = "2008" _k_eq_p3 = "2009" 	///
					_k_eq_p4 = "2010" _k_eq_p5 = "2011" _k_eq_p6 = "2012" _k_eq_p7 = "2013" 	///
					_k_eq_p8 = "2014" _k_eq_p9 = "2015"  _k_eq_p10 = "2016+") 					///
					vertical levels(95) ciopts(recast(rarea)) connect(direct) nooffsets ///
					yline(0, lcolor(black))  xline(6.4, lcolor(red)) xline(8.7, lcolor(green)) ///
					xlabel(, labsize(small) ang(45)) ylabel(-80(20)60, labsize(small) grid) ///
					ytitle("Impact on revenue per capita (in %)") name(xtevent_trend_comb_percentage, replace) ///
					legend(order(6 "Total personal tax"	- "DiD-estimate: {&beta} = `:di %3.2f ${ES_b_pcrev_total_pct}' (`:di %3.2f ${ES_se_pcrev_total_pct}')" ///
								4 "Income tax" - "DiD-estimate: {&beta} = `:di %3.2f ${ES_b_pcrev_inctax_pct}' (`:di %3.2f ${ES_se_pcrev_inctax_pct}')"  ///
								2 "Wealth tax" -  "DiD-estimate: {&beta} = `:di %3.2f ${ES_b_pcrev_weatax_pct}' (`:di %3.2f ${ES_se_pcrev_weatax_pct}')" ) ///
								row(3) ring(0) position(11) bmargin(small) size(small)) 
								
			graph export "Fig_9)-ES-total+inc+wealth-percentage_change-xtevent-pretrend.pdf", as(pdf)	replace

								* * * * *  E N D  * * * * * * 
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
